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1  Problems  Studied 

The  objectives  of  this  project  were  (1)  to  design,  implement,  and  analyze  novel  optimization 
algorithms  for  generating  the  coordinate  space  trajectories  of  molecules  that  are  undergo¬ 
ing  conformational  transformations  as  a  result  of  light-induced  excitations  and  subsequent 
relaxation  processes  and  (2)  to  improve  the  PFs  existing  Wigner-Poisson  code  to  enable 
simulations  of  devices  with  more  complex  barrier  structures. 

In  the  period  of  this  report  the  project  has  supported  three  Ph.  D.  students:  David 
Mokrauer,  who  hnished  his  Ph.  D.  in  2012,  Anne  Costolanski,  who  hnished  her  Ph.  D.  in 
2013,  and  James  Nance. 

Publications  from  the  project  include  a  book  [16],  seven  journal  papers  [5,14,17,18,27, 
29,30]  two  Ph.  D.  theses  [4,31],  two  articles  in  refereed  proceedings  [6,28],  one  submitted 
journal  paper  [35],  and  a  Sandia  technical  report  [7]. 

2  Important  Results 

2.1  Molecular  Confirmation  and  Dynamics 

The  objective  in  this  part  of  the  project  was  to  design  an  efficient  algorithm  for  simulation 
of  light-induced  molecular  transformations.  The  simulation  seeks  to  follow  the  relaxation 
path  of  a  molecule  after  excitation  by  light.  The  simulator  is  a  predictive  tool  to  see  if  light 
excitation  and  subsequent  return  to  the  unexcited  or  ground  state  will  produce  a  different 
conhguration  than  the  initial  one. 

We  simulate  the  results  of  the  excitation,  rather  than  the  excitation  itself.  The  excitation 
will  change  the  quantum  state  of  a  molecule.  Our  objective  is  to  design  software  that  will 
let  one  explore  the  possible  changes  in  a  molecule  after  a  sequence  of  excitations. 

The  goals  of  the  simulation  are  not  only  to  identify  the  end  point,  but  to  report  the 
entire  path  in  an  high-dimensional  conhguration  space  so  that  one  can  look  for  nearby  paths 
to  interesting  conhgurations  and  examine  the  energy  landscape  near  the  path  to  see  if  low 
energy  barriers  make  jumping  to  a  different  path  possible.  We  will  describe  the  hnal  version 
of  the  algorithm  from  [31].  This  version  represents  an  evolution  of  the  ideas  from  earlier 
versions  of  the  algorithm  [28-30].  The  software  has  also  been  applied  to  sensors  in  [3]. 

As  pointed  out  in  these  earlier  papers,  previous  work  on  this  type  of  problem  (see  [26,32, 
36],  for  example)  did  not  address  the  general  simulation  issue,  but  rather  were  investigations 
of  specihc  reactions  using  Gaussian  scans  of  the  surface  for  the  particular  molecule  of  interest. 
We  are  the  hrst  to  address  the  general  problem  for  more  than  a  single  degree  of  freedom. 
Our  approach  uses  the  Smolyak  sparse  interpolation  [34]  method  to  build  a  surrogate  model 


1 


of  an  expensive  molecular  dynamics  code  (in  our  particular  case  Gaussian  [11]),  and  uses 
that  model  to  drive  a  numerical  ODE  integrator. 

We  begin  in  §  2.1.1  with  a  precise  statement  of  the  problem  and  the  stages  of  model 
reduction  we  will  need  to  make  a  solution  computationally  tractable. 


2.1.1  Background 

The  task  is  to  develop  a  design  tool  which  will  simulate  molecular  changes  after  excitation. 
Ideally,  we  would  do  this  by  computing  the  ground  state,  simulating  excitation  with  a  quan¬ 
tum  chemistry  code  (in  our  case  Gaussian  [11]),  and  then  following  the  gradient  descent  path 
in  molecular  conhguration  space  (bond  lengths,  valence  angles,  and  torsion  angles). 

Molecules  of  interest  can  have  hundreds  of  atoms  and  degrees  of  freedom.  One  cannot 
vary  all  of  the  degrees  of  freedom  at  once  in  a  dynamic  simulation  because  the  simulator  is 
too  computationally  costly.  Hence,  we  must  apply  several  layers  of  model  reduction. 

The  hrst  of  these  is  to  isolate  a  few  molecular  coordinates  of  particular  interest  and  vary 
only  those  in  the  dynamic  simulation.  The  simulator  computes  an  energy  S£{p)  as  a  function 
of  a  vector  of  conhguration  variables  (torsion  and  bond  angles)  p  G  and  the  quantum 
state  £  =  0, 1, . . ..  We  split 


into  a  low-dimensional  vector  of  design  variables  x  and  the  remainder  Given  i,  we  compute 
the  energy  Ei{x)  as  a  function  of  x  alone  via 


Ee^x) 


minSf, 


(1) 


Gaussian  approximates  the  solution  of  the  optimization  problem  in  (1)  with  a  variant  of  the 
the  BEGS  [1,8,9,12,15,33]  algorithm.  The  objective  of  this  project  is  to  simulate  excitation 
of  a  molecule  from  one  quantum  state  ^  to  another  and  the  subsequent  relaxation  to  a  local 
minimum  of  E^/.  The  simulator  will  do  this  repeatedly  for  a  sequence  of  states,  beginning 
and  ending  at  £  =  0.  We  seek  a  sequence  of  states  for  which  the  initial  conhguration  at  the 
ground  state  {£  =  0)  is  diherent  from  the  one  at  the  end  of  the  sequence  of  excitations. 

For  each  £,  we  would  ideally  compute  the  local  minimum  with  the  gradient  descent 
method,  i.  e.  integrate  the  dynamics 


X  =  -VEe{x)  (2) 

with  initial  data  given  either  by  the  ground  state  conhguration  (for  £  =  0  at  the  start  of 
the  simulation)  or  by  the  result  of  excitation  from  a  local  minimum  at  another  state.  The 
evaluation  of  VEi  is  too  expensive  to  drive  the  numerical  solution  of  (2),  so  we  apply  a 
second  level  of  model  reduction  and  replace  Eg  with  a  piecewise  polynomial  interpolant. 
This  is  sufficient,  if  done  carefully,  to  solve  simple  problems  with  as  many  as  two  degrees 
of  freedom  [i.  e.  £  G  R^),  and  we  report  on  results  with  a  simple  two-dimensional  spline 
interpolation  in  [28,29].  The  problems  with  such  an  interpolation  are  that  the  number  of 
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nodes  increases  exponentially  with  the  number  of  degrees  of  freedom  (the  size  of  x)  and  one 
must  take  care  with  the  ordering  of  the  evaluation  of  the  nodes  to  make  sure  that  the  internal 
optimization  has  a  sufficiently  good  initial  iteration.  The  internal  optimization  in  Gaussian 
can  fail  if  one  does  not  pay  attention  to  the  latter  problem. 

Our  solution  to  this  complexity  problem  was  to  use  a  sparse  grid  [34]  for  the  interpolation 
nodes.  The  size  of  the  grid  grows  polynomially  with  the  dimension.  Moreover  the  grids  for 
given  orders  of  accuracy  are  nested,  enabling  us  to  perform  error  estimation  and  control. 

2.1.2  Results 

Mokrauer  defended  his  thesis  [31]  and  graduated  in  August  2012.  His  research  with  the  PI 
was  the  hrst  to  use  sparse  interpolation  as  a  surrogate  model  for  molecular  potential  energy 
functionals.  This  is  important  as  it  enables  one  to  remove  calls  to  a  quantum  chemistry 
simulator  from  a  dynamic  simulation  of  light-induced  molecular  transition  dynamics.  The 
papers  from  his  thesis  [27-30, 30]  have  all  appeared.  Mokrauer’s  LITES  code  has  been  used 
in  [3, 27, 30]  for  simulation  of  Butene,  Stilbene,  and  Azo-benzene. 

The  algorithm  in  the  LITES  code  builds  the  interpolation  adaptively,  iterpolating  over 
patches  in  conhguration  space  as  the  integration  of  the  dynamics  progresses.  While  this 
worked  well,  we  are  not  rethinking  that  approach  as  we  hnd  ways  to  interpolate  an  entire 
energy  landscape  with  a  single  interpolation. 

The  Azo-Benzene  simulation  encountered  difficulties  because  two  energy  surfaces  inter¬ 
sected.  James  Nance,  the  final  student  on  the  project,  was  working  on  making  the  evaluation 
of  the  sparse  interpolant  more  efficient  so  we  can  simultaneously  track  multiple  conforma¬ 
tional  paths  from  different  initial  points. 

When  the  project  ended,  the  PI  found  new  funding  for  Nance  and  we  hope  to  continue 
the  work.  This  will  not  only  be  useful  in  understanding  what  happens  near  local  maxima 
after  a  drop  from  an  excited  state  to  a  lower  one,  but  also  what  happens  when  a  path  crosses 
an  intersection  of  energy  surfaces.  We  must  also,  of  course,  devise  an  efficient  and  robust 
approach  to  detecting  intersections.  Our  hrst  approach  at  this  will  be  simply  to  look  for  a 
sign  change  in  the  difference  of  energies  on  the  two  paths  as  the  integration  progresses. 

Nance  and  the  PI  have  made  improvements  to  the  algorithms  in  Mokrauer’s  work.  The 
simulation  is  now  fast  enough  to  run  ensembles  of  trajectories  in  conhguration  space  to 
capture  the  ehects  of  thermal  huctuations. 

2.2  Wigner-Poisson  Solver 

In  §  2.2.1  we  will  describe  the  problem  and  the  older  work  of  the  PI  and  others.  In  §  2.2.2 
we  will  summarize  the  results  obtained  in  this  project. 

2.2.1  Background 

The  Wigner-Poisson  equation  [37]  is  an  integro-diherential  equation  that  models  the  quantum 
transport  of  electrons  in  a  semiconductor  device,  such  as  an  resonant  tunneling  diode  [2, 10, 
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38,39].  The  steady-state  equation  is  for  the  distribution  function  f{x,k),  which  depends 
on  space  x  and  momentum  k.  The  steady-state  problem  is  of  interest  in  its  own  right,  and 
for  analysis  of  dynamic  stability  of  solutions  of  the  time-dependent  problem.  We  write  the 
equation  using  the  notation  from  [38]  as 

W(f)  =  Df  +  P(f)+(^]  =0.  (3) 

\  ^  /  coll 


The  linear  term  Df  on  the  right  side  of  (3)  represents  the  kinetic  energy  effects  on  the 
distribution  and  is  given  by 

_i^dl 

27rm*  dx 

In  (4),  h  is  Planck’s  constant  and  m*  is  the  effective  mass  of  the  electron.  The  second  term, 
P{f),  is  the  nonlinear  term  in  the  equation  and  accounts  for  the  potential  energy  effects  on 
the  distribution 

4  ,  ,  , 

^  . .  (5) 


h  J-K 

The  function  T{x,  k)  is  dehned  by 


r  fix, k')T{x,k-k')dk'. 

h  J-K 


Tix,k)  = 

In  this  equation,  Uix)  is  the  electric  potential  as  a  function  of  position  and  Lc  is  the  cor¬ 
relation  length.  This  term  is  nonlinear  in  /  because  Uix)  depends  on  /  through  Poisson’s 
equation  (see  (10)  and  (12)).  The  last  term  describes  electron-electron  scattering 


\Uix  +  y)  —  Uix  —  y)]  sini2yk)dy . 


(6) 


coll 


1 

r 


f^K  fix,k')dk' 
J^Khix,k')dk' 


foix,  k)  -  fix,k) 


(7) 


In  (7),  r  is  the  relaxation  time,  and  foix,k)  is  the  equilibrium  Wigner  distribution,  /q 
is  the  solution  of  (3)  when  there  is  no  voltage  difference  (zero-bias)  across  the  device.  The 
equations  for  /o  differ  both  in  form  and  in  analytic  properties  from  those  for  the  nonzero-bias 
case. 

Boundary  conditions  are  imposed  at  the  device  edges  to  describe  the  distribution  of 
electrons  entering  the  device.  On  the  left  (x  =  0),  we  have  for  k  >  0  (electrons  with  positive 
momentum  that  are  moving  right) 


fi0,k)  =  fik) 


AnnfkBT 

Ju 


In 


^1  -|-  exp 


~  1  /  h^k‘^  y 

ksT  \8TT‘^m* 


(8) 


Similarly  on  the  right  (x  =  L)  we  specify  /  for  /c  <  0  (electrons  with  negative  momentum 
that  are  moving  left) 

fiL,k)  =  fik).  (9) 

In  (8)  and  (9),  ks  is  Boltzmann’s  constant,  T  is  the  temperature,  y  is  the  Fermi  energy  at  the 
endpoints.  The  electric  potential  Uix)  is  the  sum  of  the  potential  barrier  Ac(a;)  that  arises 
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from  the  heteroj unction  of  the  two  different  semiconductor  materials  and  the  electrostatic 
potential  u{x).  The  electrostatic  potential  is  the  solution  of  Poisson’s  equation 

Nd{x)-^  f{x,k')dk'  .  (10) 

ZTT  J-K 

In  (10),  q  is  the  charge  of  the  electron,  e  is  the  dielectric  constant,  and  Nd{x)  is  the  doping 
prohle.  Nd  is  piecewise  constant,  with  a  small  number  (<  10)  of  discontinuities.  The 
boundary  conditions  for  (10)  are 


d'^u  _ 
dx'^  e 


n(0)  =  0,u{L)  =  -Vbias,  (11) 

where  Vbias  >  0  is  the  applied  voltage  (bias).  The  potential  U  is  given  by 

U{x)  =  u{x)  +  Ac{x),  (12) 

where  Ac{x)  is  piecewise  constant  with  a  small  number  of  discontinuities. 

We  will  assume  that  the  structure  is  symmetric  about  x  =  ^,  i.  e. 

Nd{x)  =  Nd{L  -  x)  and  Ac(a;)  =  Ac(T  -  x).  (13) 

To  evalnate  the  integral  in  (6)  we  must  extend  U  outside  of  the  interval  [0,L].  We  do  this 
by  dehning  U{x)  =  U (L)  for  a;  >  L  and  U{x)  =  U (0)  for  a;  <  0.  With  this  dehnition,  U  is 
a  piecewise  smooth  function  of  x,  having  discontinuities  where  Ac  does,  /o  is  the  solution 
to  the  zero-bias  problem  (Vbias  =  0),  which  differs  from  the  nonzero-bias  (Vbias  >  0)  case  in 
that  the  collision  term  is  missing  from  the  integro-partial  differential  eqnation,  resnlting  in 

W(f)  =  Df  +  P(f)  =  0,  (14) 

and  =  0  in  the  boundary  conditions  (11)  for  the  Poisson  equation. 

Unlike  [38],  where  the  momentnm  was  nnbonnded,  older  work  [20,21,23-25,40]  of  the 
PI  and  his  collaborators  limited  k.  The  reasons  for  this  were  that  the  standard  discretiza¬ 
tion  [2,  38]  has  been  found  to  fail  to  converge  as  the  number  of  grid  points  in  momentnm 
increases  [19],  whereas  limiting  A  to  a  physically  reasonable  value  enable  not  only  nnmerical 
observations  of  convergence,  bnt  also  a  proof  [19,24]. 

2.2.2  Results 

The  work  in  the  current  project  [4-6]  began  by  correcting  some  errors  in  the  discretizations 
in  the  fortran  codes.  This  led  to  a  signihcant  improvement  in  accuracy  [5],  which  the  PI 
improved  further  with  a  grid  that  was  nonuniform  in  both  space  and  momentnm.  Costolanski 
and  the  PI  implemented  these  ideas  in  a  matlab  code.  The  matlab  code  was  promising  enough 
to  motivate  us  to  implement  the  new  discretization  in  C-I--I-  and  parallelize  the  algorithm 
using  the  Sandia  Trilinos  [13]  framework. 

Costolanski’s  hnal  C-I--I-  code  has  non-nniform  grids  in  both  space  and  momentum  [6]. 
Costolanski  gradnated  in  August  2013.  Her  thesis  reports  on  numerical  to  determine  the 
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difference  between  the  discretizations  from  [10,38]  and  those  from  [22-24],  grid  refinement 
studies,  stability  analysis,  and  continuation  results.  Her  work  corrected  several  errors  in 
previous  simulators  and,  at  least  numerically,  gave  new  insights  into  the  dynamic  stability 
of  the  Wigner-Poisson  model. 
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